One-Dimensional Stochastic Levy— Lorentz Gas. 
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We introduce a Levy-Lorentz gas in which a light parti- 
cle is scattered by static point scatterers arranged on a line. 
We investigate the case where the intervals between scatterers 
{£;} are independent random variables identically distributed 
according to the probability density function (j, (£) ~ £-( 1+ t). 
We show that under certain conditions the mean square dis- 
placement of the particle obeys (x 2 (t)) > Ct 3 ~ J for 1 < 
7 < 2. This behavior is compatible with a renewal Levy 
walk scheme. We discuss the importance of rare events in the 
proper characterization of the diffusion process. 



PACS numbers: 02.50.-r, 05.40. +j, 05.60,+w 



I. INTRODUCTION 

In recent years there has been a growing interest in 
anomalous diffusion defined by 



(x 2 ) = D s t s 



(1) 



and S > 1 [Q-Q]. Such a behavior was found in chaotic 
diffusion in low dimensional systems |5|j^| , tracer diffusion 
in a rotating flow |t]] , N body Hamiltonian dynamics || , 
Lorentz gas with infinite horizon |l],|9| and diffusion in egg 
crate potentials |T(|. In all these examples one observes 
long ballistic nights in which the diffusing particle moves 
at a constant velocity. The transport is characterized by 
a distribution of free flight times which follows a power 
law decay. These processes have been usually analyzed 
using the Levy walk framework (see more details below) 

ill- „ 

It has been recently suggested by Levitz |0] that three- 
dimensional molecular Knudsen diffusion, at very low 
pressures, inside porous media can be described by Levy 
walks. It has been also shown that pore chord distri- 
butions, in certain three-dimensional porous media de- 
cay as a power law, at least for several length scales. 
Hence one can anticipate that a light test particle in- 
jected into such a medium may exhibit a Levy walk. This 
has motivated the investigation of a fractal Lorentz gas. 
Levitz [|ll| has simulated trajectories of a light particle 



reflected from a three-dimensional intersection of a four- 
dimensional Weierstrass-Mandelbrot hyper surface, and 
found an enhanced Levy type diffusion. 

Here we investigate a one dimensional stochastic 
Lorentz gas which we call Levy-Lorentz gas. In this 
model a light particle is scattered by a fixed array of iden- 
tical scatterers arranged randomly on a line. Upon each 
collision event the light particle can be transmitted (or 
reflected) with probability T (or R = 1 — T). We inves- 
tigate the case when the intervals between the scatterers 
are independent identically distributed random variables 
with a diverging variance. 

We find: (a) a lower bound for the mean square dis- 
placement which is compatible with the Levy walk model, 
and (b) that the generalized diffusion coefficient D$ is 
very sensitive to the way the system has been prepared 
at time t — 0. As expected, we show that the transport 
is not Gaussian. In systems that exhibit normal diffu- 
sion, the contribution from ballistic motion, x 2 — v 2 t 2 , 
is important only for short times; here we show that the 
ballistic motion cannot be neglected even at t — > oo. The 
ballistic paths contribute to the generalized diffusion co- 
efficient D$ exhibiting a behavior different than normal. 



II. MODEL AND NUMERICAL PROCEDURE 

Assume a light particle which moves with a constant 
speed (v = ±1) among identical point scatterers arranged 
randomly on a line. Upon each collision, the probabil- 
ity that the light particle is transmitted (reflected) is T 
(R = 1 — T). The intervals between scattering points, 
fi > with (« = •••, — n, • • • , — 1, 0, 1, • • •), are indepen- 
dent identically distributed random variables described 
by a probability density function /i(£). An important 
random variable is Xf defined to be the distance between 
the initial location of the light particle (x = 0) and the 
first scatterer in the sequence located at x > 0. The ran- 
dom variable Xf is described by the probability density 
function h(xf). A set of scatterers (black dots) is given 
schematically by: 
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where the open circle represents the light test particle 
at time t = 0. We consider the case when for large £, 
M(£) ~ <?T (1+7 \ with < 7 < 2. Thus the variance of 
the length intervals diverges. A realization of the 
scatterers is shown in Fig. [I], for the case 7 = 3/2. We 
observe large gaps which are of the order of the length of 
the system. 

The case for which the variance converges has been in- 
vestigated thoroughly in [^3^171, resulting in: (i) a nor- 
mal Gaussian diffusion as expected from the central limit 
theorem, and (ii) a 3/2 power law decay in t of the ve- 
locity autocorrelation function. 

Along this work we present numerical results for the 
case 7 = 3/2 and T — 1/2. We use the following nu- 
merical procedure. First we generate a set of scatterers 
on a one dimensional lattice with a lattice spacing equal 
unity. Using a discrete time and space iteration scheme 
we find an exact expression for the probability of finding 
the particle on x at time t, p(x, t\x = 0, t = 0), given that 
at t = the particle is on x = 0. The initial location of 
the particle is determined using equilibrium initial con- 
ditions (see details below). The initial velocity is v = 1 
or v — — 1 with equal probabilities. p{x,t\x = 0, t = 0) 
depends on the realization of disorder we have generated 
in the first step. This procedure is repeated many times. 

In appendix A we explain how we generated the ran- 
dom intervals When 7 = 3/2 the mean (£) = 
Ji} 00 ls nm te while the second moment (£ 2 ) = 00. 
Since \v\ = 1 the characteristic microscopic time scale is 
(£) which is referred to as the mean collision time, and 
our simulations are for times t ~ 1000 (£). For our choice 
of parameters (£) ~ 4 (see more details in Appendix A). 

III. RESULTS 

Let us analyze our one-dimensional Levy-Lorentz 
model using the Levy walk approach |2[-|i|JqJlq| . Levy 
walks describe random walks which exhibit enhanced dif- 
fusion and are based on the generalized central limit theo- 
rem and Levy stable distributions |l^| . Briefly, a particle 
moves with a constant velocity » = +1 or 11 = -1 and 
then at a random time t\ its velocity is changed. Then 
the process is renewed. Each collision is independent 
of the previous collisions. The times between collision 
events {r^} are assumed to be independent identically 
distributed random variables, given in terms of a proba- 
bility density function q (r) . One might expect that the 
dynamics of Levy-Lorentz gas can be analyzed using the 
Levy walk renewal approach with q(r) ~ t~^ 1+1 \ for 
large r and < 7 < 2, which leads to 



( t 3 ~^, 1 < 7 < 2 
(x 2 ) ~ { (2) 
I t 2 , 0<7<1. 

For 7 > 2 one finds normal diffusion. It is clear that 
the renewal Levy walk approach and the Levy-Lorentz 
gas are very different. Within the Levy-Lorentz gas colli- 
sions are not independent and correlations are important. 
Hence it is interesting to check whether the renewal Levy 
walk model is suitable for the description of the Levy- 
Lorentz gas. In this context it is interesting to recall that 
Sokolov et al (2(J have shown that correlations between 
jumps in a Levy flight in a chemical space destroy the 
Levy statistics of the walk. 
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FIG. 1. A realization of a set of scatterers with 7 = 3/2 
exhibiting gaps on many scales. The horizontal axis is the 
x coordinate. All along this work we consider dimensionless 
units. 

We consider a continuum model to derive our ana- 
lytical results; the generalization to the lattice case is 
straightforward. Let (p (x, t\x = 0, t = 0))dx be the prob- 
ability, averaged over disorder, of finding the test parti- 
cle at time t, in the interval (x, x + dx). Initially, at time 
t = 0, the particle is at x — 0, and there is an equal prob- 
ability of the particle having a velocity v = +lovv = — 1. 
Figs. ^| and ^present numerical simulations which show 
(p (x, t\x = 0, t = 0)). One can see that in addition to the 
central peak on x = 0, two other peaks appear at loca- 
tions x = ±t. These peaks, known as ballistic peaks, were 
observed in a similar context in other systems exhibiting 
enhanced diffusion P, pdpl| ] . The peaks are stable on the 
time scale of the numerical simulation. The height of 
these peaks decays with time, and according to our finite 
time numerics the central peak and the ballistic peaks 
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decay according to the same power law when 7 = 3/2. 

Let us analyze analytically the time dependence of the 
ballistic peaks and calculate their contribution to the 
mean square displacement. Since in our model \v\ = 1 it 
is clear that 



(p(x,t\x = 0,t = 0)) =0, for|x|>t. 



(3) 



We decompose the ensemble averaged probability density 
into two terms 

(ja(x,t\x = 0,t = 0)) = 

(p (x, t\x = 0, t = 0)> + ^Q b (t) [S(x + t) + S(x~ t)] . 

(4) 

The first term on the RHS, (p(x,t\x = 0,t= 0)), is the 
probability density of finding the light particle at |x| < t. 
Qb if) is the probability of finding the light particle at 
x = t (x = —t) if initially at x = and its velocity is +1 
(— 1). The left-right symmetry in Eq. (Q) means that we 
have used the symmetric initial condition (i.e., v = +1 
or v = — 1 with equal probabilities) and the assumption 
that the system of scatterers is isotropic in an averaged 
sense. Using a similar notation we write 



FIG. 2. A histogram presenting the (p(x, t\x = 0,t = 0)) 
versus x for time t = 1024, 7 = 3/2 and T = 1/2. Notice the 
ballistic peaks of the propagator at x = it. The average is 
over 1.8* 10 5 realizations of disorder. The bin length is unity. 



(x 2 ) 



(5) 



where (x 2 )b is the ballistic contribution to the mean 
square displacement. From Eq. (H) we have 



h(t)t 2 <(x 2 {t))<t 2 



(6) 



The upper bound is an obvious consequence of the fact 
that \v\ = 1. The lower bound, found using (x 2 (t))b < 
(x 2 (t)}, is of no use when all moments of \i (£) converge, 
since then Qb (t) usually decays exponentially for long 
times. Eq. (Q) is useful when the moments of di- 
verge, a case we consider here. 
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FIG. 3. The same as Fig. 

^| for time t = 4096. The probability of finding a ballistic 
path, (p(x = ±t,t\x — 0,t — 0)} ~ 0.004, is small but yet of 
statistical significance when (x 2 (t)) is calculated. 

To find Qb (t) consider a test particle which is initially 
of velocity +1 and located at x = 0. The probability it 
reaches x = t, at time t, is T r where r is the number of 
scatterers in the interval of length (0,t). Hence, 
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Q b (t)=Y J T r G r (t), 



(7) 



and G r (t) is the probability of finding r scatterers in 
(0, t). G r (t) can be calculated in terms of and of 
h(xf). In appendix B we use renewal theory to calculate 
the Laplace t — > u transform of Q b (t) 



W s 1 (T-l)h(u) 



[1 - Tfx (u)] ■ 



(8) 



When T = 1, Qb (u) = l/u, as expected from a trans- 
mitting set of scatterers. In deriving Eq. (||) we have 
used the model assumptions that the intervals {^} are 
statistically independent and identically distributed. 

The function h(xf) depends on the way the system of 
scatterers and light particle are initially prepared. Con- 
sider the following preparation process. A scatterer is 
assigned at the location x = —L (eventually L — > oo), 
then random independent length intervals are generated 
using the probability density These length intervals 
determine the location of scatterers on the line. When 
the sum of the length intervals exceeds 2L the process 
is stopped. As mentioned, at time t = the light par- 
ticle is assigned to the point x = 0. When the mean 
distance between scatterers (£) = J °° £^(£)^£ converges 
(i.e., 1 < 7), and L — > 00 then according to fl9|j22|l 



h{x f ) 



(0 



(9) 



which is standard in the context of the Lorentz gas when 
the moments of /u(£) converge |L3|,[L5). This type of initial 
condition is called equilibrium initial condition. When 
1 < 7 < 2, Eq. (||) implies that h(xf) ~ (£/)~ 7 an d 
hence (xf) — x fh(xf)dx / — > 00. At first sight this 
divergence might seem to be paradoxical, since the mean 
distance between scatterers, (£), converges. We notice 
however that the point x = has a higher probability 
to be situated in a large gap. Hence, statistically the 
interval £0 is much larger than the others and in our case 
(xj) = 00. Eq. (|^) implies that on average one has to 
wait an infinite time for the first collision event. 
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FIG. 4. The probability to find the first scatterer at a dis- 
tance Xf from the origin. Here the average is over 3 * 10 
realizations, and half the length of the system is L — 10 . We 
use a bin of length 32 (dimensionless units). The solid curve 
is the theoretical prediction, Eq. (^|), with no fitting param- 



eters. For large 27, h(x/) 
that (xf) diverges. 



x, 1 and 7 = 3/2, which implies 



In numerical simulations the system's length L is finite, 
so that Eq. (^) is only an approximation which we expect 
to be valid for Xf « L. However, if we observe a system 
for time t « L the boundary condition is not expected 
to influence the anomalous dynamics. We have generated 
numerically many random systems, using /i(£) ~ £-(!+7) 
and 7 = 3/2. As shown in Fig. |[ h(xf) ~ (x/)~ 7 as 
predicted in Eq. (||). 

We consider the small u expansion, of the Laplace 
transform of /i(£), 



(10) 



where 1 < 7 < 2 and a is a constant. Using a Tauberian 
theorem and Eqs. (0), (^|) it can be shown that for long 
times t 



r(2- 7 ) \(0 



Q b (t) = 

2a (7-I) T 



F(2- 7 ) l-T V(0 



Inserting Eq. (|Tl|) in Eq. @ we find 



f 



3-7 



777 <{x'{t))<t' 



(11) 



(12) 



r(2- 7 ) V(0, 

This bound demonstrates that the diffusion is enhanced, 
namely the mean square displacement increases faster 
than linearly with time. 
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In Fig. |[ we present the mean square displacement of 
the light particle obtained by numerical simulations for 
the case 7 = 3/2. We see that for the chosen values of 
parameters the asymptotic i 3-7 behavior can be observed 
for times which are accessible on our computer. Fig. |B| 
clearly shows that the ballistic contribution (x 2 )h to the 
mean square displacement (x 2 ) is significant. Notice that 
our numerical results are presented for times which are 
much larger than the mean collision time (£) ~ 4. 

In Fig. ^| we show the probability of finding a ballistic 
path, namely, the probability of finding the light particle 
at time t at x = +t, or at x = —t. By definition these 
probabilities are equal to Qb(t)/2. We observe the i 1-7 
behavior of Qb(t), Eq. (|ll"l), with which our lower bound 
was found. The fact that the probability of finding the 
particle at x = t is equal to the probability of finding 
the particle at x = — t means that our system is isotropic 
in an averaged sense. This is achieved by choosing large 
values of L. 



Only under this condition the behavior in Eq. (|12j) is 
expected to be valid. 

The lower bound in Eq. (|l^) is compatible with the 
renewal Levy walk approach Eq. (||). Other stochastic 
models |l]j23|] for enhanced diffusion based on Levy scal- 
ing arguments predict 



(x 2 ) - t 2h for 1< 7 < 2 



(14) 



which is different from Eq. (g). This approach is based 
upon a fractional Fokker-Planck equation (FFPE) 



dp (x, t) 
dt 



= L» 7 V>(x,i) 



(15) 



used in p3fl to predict an enhanced diffusion. The non- 
local fractional operator in Eq. ( |l5| ) is defined in Fourier 
k space according to the transformation V 7 — > — |/c| 7 . 
Our findings here show that Eq. (14) does not describe 
the dynamics of the Levy-Lorentz gas, since 3 — 7 > 2/7 
for 1 < 7 < 2. 




Log M (t) 



100 1000 

L°e, (t) 



FIG. 5. log 10 [(z 2 )] versus log 10 (i). The points are numer- 
ical results. The straight curve is the asymptotic behavior of 
© (i.e., (x 2 ) b ). We use 7 = 3/2 and 



the lower bound, Eq. 

so (x 2 ) > Ct i/2 . 



The lower bound in Eq. ( |l2] ) does not depend on the 
transmission coefficient T. Thus, even when all the scat- 
terers are perfect reflectors, with R = 1, the diffusion is 
enhanced. Large gaps which are of the order of the length 
t are responsible for this behavior. The transmission co- 
efficient has an important role in determining what is the 
asymptotic time of the problem. The condition that the 
first term in Eq. (Jill) dominates over the second reads: 



T 



(13) 



FIG. 6. The probability to find the light particle at time 
t at x = +t (stars) and at x = — t (dots) versus time. The 
solid curve is the theoretical prediction, Eq. ([y]) (no fitting 
parameters) which gives Qb{t)/2 ~ t' 1 " 7 ', with 7 = 3/2. 

Consider now the case when the light particle is ini- 
tially located at a scattering point. Such an initial con- 
dition is called non equilibrium initial condition. Under 
this condition h(xf) = fJ>{xf) instead of Eq. (Q). This 
means that the particle has to wait an average time (£) 
before the first collision event instead of the infinite time 
when the equilibrium initial conditions were used. Using 
Eqs. (§), (|) and @, and fi(u) = 1 - {An) 1 + ■■■ for 
< 7 < 1 and small u, we find 



5 



(x 2 )> 



2-7 



1-T T(2-7 



1 < 7 < 2 
< 7 < 1. 



(16) 



For 1 < 7 < 2 the bound differs from the i 3 ~ 7 found in 
Eq. (|l2|), where we chose h(xf) according to Eq. (£). 
Thus the ballistic contribution (x 2 )b, defined in Eq. (5), 
behaves differently for the two ensembles even when t — > 
oo. This is very different from regular Lorentz gases, 
which in the limit t — > oo are not sensitive to the choice 
of h(x f ). 

Finally, Fig. [?] shows the behavior of the correla- 
tion function (p(x = 0, t\x = 0, t = 0)) obtained from 
the numerical simulation with equilibrium initial condi- 
tions. We observe a i -1 / 2 decay of the correlation func- 
tion. This behavior is compatible with standard Gaus- 
sian diffusion which gives the well known t~ d l 2 result 
in d dimensions. We find this behavior for time scales 
which arc much larger than the mean collision time (£), 
however we have no proof that this behavior is asymp- 
totic. On the other hand the Levv_ walk model predicts 
(p(x = 0, t\x = 0, t = 0)) ~ * 



gas is compatible with the Levy walk framework and not 
with the FFPE. 

(b) Ballistic contributions to the mean square displace- 
ment are important even for large times. 

(c) The ballistic peaks at x = +t and x = —t can be 
analyzed analytically. They decay as power laws. 

(d) The way in which the system is prepared at t = (i.e., 
equilibrium versus non equilibrium initial conditions) de- 
termines the behavior of the ballistic peaks. Since these 
peaks contribute to the mean square displacement even 
at large times , we conclude that the diffusion coefficient 
Dg is sensitive to the way the system is prepared. 

In our work we considered an initial condition v = 1 
or v — — 1 with equal probabilities. It is clear that if we 
assign a velocity v — +1 to the light particle, at t = 0, 
(j)(x,t\x — 0,t — 0)) will never become symmetric, even 
approximately. Instead of the three peaks in Fig. || one 
will observe only two peaks one at x — and the other 
at x = +t. 

The reason for these behaviors in the Levy-Lorentz gas 
stems from the statistical importance of ballistic paths. 
This is different from the systems in which diffusion is 
normal in which these paths are of no significance at long 
times. Thus, similarly to Newtonian dynamics, the sys- 
tem exhibits a strong sensitivity to initial conditions. 

Experiments measuring diffusion phenomena usually 
sample data only in a scaling regime (e.g., < 
x < \[U\i ). Rare events where the diffusing particle 
is found outside this regime are many times assumed to 
be of no statistical importance. Here we showed that for 
the Levy-Lorentz gas rare events found in the outer most 
part of (p(x, t\x = 0, t = 0)) are of statistical importance. 



Note added in proof. Recently related theoretical work 
on enhanced diffusion was published ^3 

Acknowledgment We thank A. Aharony, P. Levitz, 
I. Sokolov and R. Metzler for helpful discussions. 
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FIG. 7. log 10 [{p(x = Q,t\x = 0,t = 0))] versus log 10 (i). 
The points are numerical results. The straight curve is a 
fit exhibiting the t" 1 ^ 2 behavior. 



IV. SUMMARY AND DISCUSSION 

In this work we have considered a one dimensional 
Levy-Lorentz gas. We have shown that: 
(a) the mean square displacement in the Levy-Lorentz 



A. Appendix A 

As mentioned we use a lattice model for the simulation 
so that £ is an integer. We use the transformation 

tan I 



£ = INT 



1. 



(17) 



Here I = INTjzj is the integer closest to z satisfying 
I < z. In Eq. ( |17| ) u is a random variable distributed 
uniformly according to 



< Un 



< u < u n 



< 1, 



(18) 



where u rn i n and u max are cutoffs. It is easy to generate 
the random variable u on a computer. The probability 
to find an interval of length £ is, 



G 



MO 



Hc{y)dy 



(19) 



with 



Hc(y) = < 



y < y mm 

7rA l+y 2 ~i ymin < J/ < 2/max 

(20) 

2/ > 



Here 



tan 



1/7 



y-max — 



tan 



l/ 7 

(21) 



are the cutoffs of fi c (y). When w m i„ = and u max = 1 
we have y mm = and y max = oo. In Eq. @ 
A = « mal — u m i„ determines the normalization condi- 
tio 11 JrT ^{y)dy = 1. 

To derive Eqs. (|l9|)-(^l|) we use the transformation y = 

[tan (wr/2)] 1 / 7 , and then /u c (y) = g (u) |du/<iy|, where 
g (u) is the uniform probability density of it. 
For large £ we find 



m(0 



2x^-1-7 C < t 



(22) 







with £ rocM = INT{[tan(u mQ2; 7r/2)] 1/7 } + 1. When 

u mra = 1 the second moment of diverges. 

In our numerical simulations we consider u m i„ = 1/2, 
ii mal = 1 and 7 = 3/2. Then (£) ~ 4.031 and for large 
£, we find 



m(0 ~ (6A)r 5/2 e » i. 



B. Appendix B 



(23) 



The calculation of Qb(u) can be found in |T^.^|. The 
probability that the interval (0, t) is empty is 



G (t) = 1 - / h(r)dT 
Jo 



(24) 
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and in Laplace space Go (it) = [1 — h(u)]/u. The Laplace 
transform of G r (t) for r > 1 is found using convolution 



G r (u) = /i(u)A r_1 (u)W(M) 



(25) 



W(t) = 1 — J * fj,(£)d£ is the probability that an interval 
of length (0, t) is empty, given that a scatterer occupies 
0~. In Laplace space W(u) = [1 — jl(u)]/u. Using Eqs. 
©, © and (0) we find (§. 
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